Load packages

library(tidyverse)
library(gtsummary)
library(flextable)
library(corrplot)
library(pacman)
library(devtools)
library(qgcomp)
library(forcats)
library(cowplot)
library(broom)
library(bkmr)
library(qgcompint) 
library(kohonen)
library(sommix)
library(prettydoc)

pacman::p_load(future, future.apply, reshape2, MASS, gWQS)

Descriptive Statistics

Data Exploration

Table 1

Distribution of sociodemographic characteristics and Child Behavior Checklist (CBCL) outcomes among Black mother-child pairs in metropolitan Atlanta, Georgia.

# Select variables
table_1 <- mix %>% dplyr::select(`Maternal Age`, `Maternal Education`, `Marital Status`, Parity,`Pre-pregnancy BMI`,`Substance Use`, `Insurance Type`, `Income to Poverty Ratio`, `Infant Sex`, `Mean Child Age`, Internalizing, Externalizing) 

# Set theme
theme_gtsummary_journal(journal = "lancet")

# Table
table_1 <- table_1 %>% 
  tbl_summary(type = all_dichotomous() ~ "categorical", statistic = all_continuous() ~ "{mean} ({sd})", digits = all_continuous() ~ 2) %>%
  modify_header(label = "**Demographic**")  %>%
  bold_labels()

table_1
Demographic N = 1591
Maternal Age 24·97 (4·46)
Maternal Education
    Less than some college 83 (52%)
    College degree or higher 76 (48%)
Marital Status
    Never married 75 (47%)
    Married 84 (53%)
Parity
    No prior births 65 (41%)
    1 or more prior births 94 (59%)
Pre-pregnancy BMI
    Underweight or normal 64 (40%)
    Overweight 29 (18%)
    Obese 66 (42%)
Substance Use
    No 90 (57%)
    Yes 69 (43%)
Insurance Type
    Low-income Medicaid 69 (43%)
    Pregnancy Medicaid 58 (36%)
    Private 32 (20%)
Income to Poverty Ratio
    <100% 76 (48%)
    100% - 199% 52 (33%)
    >200% 31 (19%)
Infant Sex
    Male 81 (51%)
    Female 78 (49%)
Mean Child Age 2·82 (0·64)
Internalizing 7·10 (5·29)
Externalizing 11·37 (7·62)
1 Mean (SD); n (%)

Table S1

Distributions of per- and polyfluoroalkyl substances (PFAS, ng/mL) and polybrominated diphenyl ethers (PBDEs, pg/mL), measured in serum samples obtained during the first trimester, among Black pregnant women in metropolitan Atlanta, Georgia.

# Select variables
table_s1 <- mix %>% dplyr::select(PFHXS, PFOS, PFOA, PFNA, PBDE47, PBDE99, PBDE100)
table_s1_PFAS <- atl_chemicals %>% dplyr::select(PFHXS, PFOS, PFOA, PFNA) %>% drop_na(PFHXS, PFOS, PFOA, PFNA)
table_s1_PBDE <- atl_chemicals %>% dplyr::select(PBDE47, PBDE99, PBDE100) %>% drop_na(PBDE47, PBDE99, PBDE100)

# Build function 
dist.func <- function(x) {
  c(
    geomean = round(exp(mean(log(x), na.rm = TRUE)), 2),
    geosd = round(exp(sd(log(x), na.rm = TRUE)), 2),
    min = round(min(x, na.rm = TRUE), 2),
    max = round(max(x, na.rm = TRUE), 2),
    median = round(median(x, na.rm = TRUE), 2)
  )
}

# Apply 
table_s1 <- as.data.frame(do.call(rbind, lapply(table_s1, dist.func)))
table_s1 <- table_s1 %>%
  mutate(sample = "Analytic")
table_s1_PFAS <- as.data.frame(do.call(rbind, lapply(table_s1_PFAS, dist.func)))
table_s1_PFAS <- table_s1_PFAS %>%
  mutate(sample = "Full")
table_s1_PBDE <- as.data.frame(do.call(rbind, lapply(table_s1_PBDE, dist.func)))
table_s1_PBDE <- table_s1_PBDE %>%
  mutate(sample = "Full")

table_s1 <- rbind(table_s1_PFAS, table_s1, table_s1_PBDE)
# Table
table_s1 <- table_s1 %>% 
mutate(Chemical = rownames(table_s1)) %>%
mutate("Geometric Mean (SD)" = paste0(geomean, " (", geosd, ")")) %>%
dplyr::select(Chemical, sample, `Geometric Mean (SD)`, min, median, max) %>%
  as_grouped_data("sample") %>%
  flextable() %>%
  set_header_labels(
  geomean = "Geometric Mean",
  geosd = "Geometric SD",
  min = "Minimum",
  max = "Maximum",
  median = "Median") %>% 
  font(fontname = "Tahoma") %>%
  autofit() %>%
  theme_alafoli()

table_s1

sample

Chemical

Geometric Mean (SD)

Minimum

Median

Maximum

Full

PFHXS

1.16 (2.06)

0.08

1.24

6.17

PFOS

1.89 (2.5)

0.01

2.12

12.42

PFOA

0.62 (2.36)

0.01

0.70

4.42

PFNA

0.26 (2.35)

0.01

0.30

2.27

Analytic

PFHXS1

1.07 (2.07)

0.14

1.27

4.50

PFOS1

2.37 (1.97)

0.14

2.53

12.42

PFOA1

0.8 (2.12)

0.03

0.89

4.42

PFNA1

0.33 (2.04)

0.01

0.34

2.27

PBDE47

92.7 (1.98)

20.65

87.48

441.86

PBDE99

23.07 (2.21)

7.81

23.26

308.29

PBDE100

14.6 (2.85)

2.81

17.28

111.47

Full

PBDE471

89.84 (2.08)

15.12

83.86

1,382.08

PBDE991

23.1 (2.35)

7.81

22.11

341.98

PBDE1001

13.32 (2.99)

2.81

15.44

289.20

Figure S1

Pearson correlation for the association between natural log transformed per- and polyfluoroalkyl substances (PFAS, ng/mL) and polybrominated diphenyl ethers (PBDEs, pg/mL), measured in first trimester serum samples, and internalizing and externalizing behaviors averaged across ages 1-5 years among Black mother-child pairs in metropolitan Atlanta, Georgia.

# Correlation matrix and plot
mix_cor <- mix %>% dplyr::select(PFHXS, PFOS, PFOA, PFNA, PBDE47, PBDE99, PBDE100, Internalizing, Externalizing)

cor <- cor(mix_cor, method = "pearson")

figure_s3 <- corrplot(cor, method = "color", tl.col = "grey25", tl.cex = 6, cl.cex = 6, mar=c(0,0,1,0)) %>%
corrRect(c(1,5,8, 9))

Regression

Multiple Linear Regression

Model

# Total Sample
# Assign predictors, outcomes
predictors <- c("PFHXS.log", "PFOS.log", "PFOA.log", "PFNA.log", "PBDE47.log", "PBDE99.log", "PBDE100.log")
outcomes <- c("Internalizing", "Externalizing")

# Loop regression
models <- list()
for (i in 1:length(outcomes)) {
  for (j in 1:length(predictors)) {
    model_name <- paste0(outcomes[i], " ~ ", predictors[j])
    formula <- as.formula(paste(outcomes[i], " ~ ", predictors[j], " + `Mean Child Age`  + `Infant Sex` + `Marital Status` + `Pre-pregnancy BMI` + `Maternal Age` + `Maternal Education` + Parity + `Substance Use`"))
    model <- lm(formula = formula, data = mix)
    models[[model_name]] <- tidy(model, conf.int = TRUE) %>%
      mutate(outcome = outcomes[i])
  }
}
chemicals_adj <- do.call(rbind, models) 

# Format 
reg <- chemicals_adj %>% 
  filter(term %in% c("PFHXS.log", "PFOS.log", "PFOA.log", "PFNA.log", "PBDE47.log", "PBDE99.log", "PBDE100.log"))  %>%
  dplyr::select(outcome, term, estimate, conf.low, conf.high) %>%
  mutate(sex = "All")

### By sex

# Male
## Filter for males
mix_m <- mix %>% filter(`Infant Sex` == "Male")

models <- list()
for (i in 1:length(outcomes)) {
  for (j in 1:length(predictors)) {
    model_name <- paste0(outcomes[i], " ~ ", predictors[j])
    formula <- as.formula(paste(outcomes[i], " ~ ", predictors[j], " + `Mean Child Age` + `Marital Status` + `Pre-pregnancy BMI` + `Maternal Age` + `Maternal Education` + Parity + `Substance Use`"))
    model <- lm(formula = formula, data = mix_m)
    models[[model_name]] <- tidy(model, conf.int = TRUE) %>%
      mutate(outcome = outcomes[i])
  }
}
chemicals_adj <- do.call(rbind, models) 

reg_male <- chemicals_adj %>% 
  filter(term %in% c("PFHXS.log", "PFOS.log", "PFOA.log", "PFNA.log", "PBDE47.log", "PBDE99.log", "PBDE100.log"))  %>%
  dplyr::select(outcome, term, estimate, conf.low, conf.high) %>%
  mutate(sex = "Male")

# Female
## Filter for females
mix_f <- mix %>% filter(`Infant Sex` == "Female")

models <- list()
for (i in 1:length(outcomes)) {
  for (j in 1:length(predictors)) {
    model_name <- paste0(outcomes[i], " ~ ", predictors[j])
    formula <- as.formula(paste(outcomes[i], " ~ ", predictors[j], " + `Mean Child Age` + `Marital Status` + `Pre-pregnancy BMI` + `Maternal Age` + `Maternal Education` + Parity + `Substance Use`"))
    model <- lm(formula = formula, data = mix_f)
    models[[model_name]] <- tidy(model, conf.int = TRUE) %>%
      mutate(outcome = outcomes[i])
  }
}
chemicals_adj <- do.call(rbind, models) 

reg_female <- chemicals_adj %>% 
  filter(term %in% c("PFHXS.log", "PFOS.log", "PFOA.log", "PFNA.log", "PBDE47.log", "PBDE99.log", "PBDE100.log"))  %>%
  dplyr::select(outcome, term, estimate, conf.low, conf.high) %>%
  mutate(sex = "Female")

Figure 1

Associations between natural log transformed per- and polyfluoroalkyl substances (PFAS, ng/mL) and polybrominated diphenyl ethers (PBDEs, pg/mL), measured in serum samples obtained during the first trimester, and internalizing and externalizing behaviors averaged across ages 1-5 years among Black mother-child pairs in metropolitan Atlanta, Georgia. PFAS in ng.mL and PBDEs in pg/mL.

Models are adjusted for infant sex, mean child age, maternal age, maternal education, marital status, early pregnancy BMI, parity, and substance use

# Format
reg$outcome <- factor(reg$outcome, levels = c("Internalizing", "Externalizing"))
reg$term <- factor(reg$term, levels = c("PFHXS.log", "PFOS.log", "PFOA.log", "PFNA.log", "PBDE47.log", "PBDE99.log", "PBDE100.log"))
chemical_label <- c("PFHXS", "PFOS", "PFOA", "PFNA", "BDE-47", "BDE-99", "BDE-100")
names(chemical_label) <- c("PFHXS.log", "PFOS.log","PFOA.log", "PFNA.log", "PBDE47.log", "PBDE99.log", "PBDE100.log")

# Figure
figure_1 <- ggplot(reg, aes(x = estimate, y = term, color = outcome)) +
  geom_pointrange(aes(xmin = conf.low, xmax = conf.high), size = 4, linewidth = 4) +
  geom_vline(xintercept = 0, color = "grey25", linetype = "dashed", alpha = 0.5) +  
  facet_grid(~outcome, scales = "free") +
  xlab("Beta (95% Confidence Interval)") +
  ylab("Chemical") +
  scale_x_continuous(n.breaks=4) +
  scale_y_discrete(labels = chemical_label) +
  scale_color_manual(values = c("grey25", "grey55")) +
  coord_flip() +
  theme_classic(base_family = "AppleGothic") +
  theme(title = element_text(size = 50, color = "grey25"),
        axis.title.x = element_text(size = 50, color = "grey25"),
        axis.title.y = element_text(size = 50, color = "grey25"),
        axis.text.x = element_text(size = 50, color = "grey25", angle = 45, vjust = 0.6),
        axis.text.y = element_text(size = 50, color = "grey25"),
        strip.background = element_rect(fill="white"),
        strip.text = element_text(colour = 'grey25', face = "bold", size = 50), 
        legend.position = "none", 
        strip.background.x = element_blank(),
        panel.grid.major = element_blank(), panel.grid.minor = element_blank())

figure_1

Table S2

Linear regression associations (overall, child sex specific) between natural log transformed per- and polyfluoroalkyl substances (PFAS, ng/mL) and polybrominated diphenyl ethers (PBDEs, pg/mL), measured in serum samples obtained during the first trimester, and internalizing and externalizing behaviors averaged across ages 1-5 years among Black mother-child pairs in metropolitan Atlanta, Georgia.

Total sample models are adjusted for infant sex, mean child age, maternal age, maternal education, marital status, early pregnancy BMI, parity, and substance use. Sex specific models are adjusted for mean child age, maternal age, maternal education, marital status, early pregnancy BMI, parity, and substance use.

# Combine regression results
reg <- rbind(reg, reg_female, reg_male)

# Table
table_s2 <- reg %>%
  mutate_if(is.numeric, round, digits = 2) %>%
  mutate(Beta_CI = paste0(estimate, " (", conf.low, ", ", conf.high, ")"))  %>%
 mutate(term = recode(term, 
                     "PFHXS.log" = "PFHXS",
                     "PFOA.log" = "PFOA",
                     "PFNA.log" = "PFNA",
                     "PFOS.log" = "PFOS",
                     "PBDE47.log" = "BDE-47",
                     "PBDE99.log" = "BDE-99",
                     "PBDE100.log" = "BDE-100")) %>%
  dplyr::select(term, outcome, sex, Beta_CI) %>%
  pivot_wider(names_from = outcome, values_from = (Beta_CI)) %>%
  group_by(term, sex) %>%
  summarise(
    Internalizing = Internalizing,
    Externalizing = Externalizing) %>%
  as_grouped_data(groups = "term") %>%
  flextable() %>%
  add_body_row(values = c("","" ,"Beta (95% CI)", "Beta (95% CI)")) %>%
  font(fontname = "Tahoma") %>%
  set_header_labels(term = "Chemical", sex = "Sex") %>%
  autofit() %>%
  theme_alafoli()

table_s2

Chemical

Sex

Internalizing

Externalizing

Beta (95% CI)

Beta (95% CI)

PFHXS

All

-0.17 (-0.33, 0)

-0.21 (-0.36, -0.05)

Female

-0.12 (-0.35, 0.11)

-0.17 (-0.37, 0.02)

Male

-0.17 (-0.42, 0.07)

-0.22 (-0.48, 0.03)

PFOS

All

-0.1 (-0.26, 0.06)

-0.21 (-0.36, -0.05)

Female

-0.01 (-0.25, 0.23)

-0.13 (-0.33, 0.06)

Male

-0.16 (-0.39, 0.06)

-0.24 (-0.48, -0.01)

PFOA

All

-0.12 (-0.28, 0.04)

-0.17 (-0.33, -0.01)

Female

-0.02 (-0.36, 0.32)

-0.13 (-0.42, 0.15)

Male

-0.18 (-0.39, 0.02)

-0.21 (-0.43, 0)

PFNA

All

-0.15 (-0.32, 0.01)

-0.21 (-0.36, -0.05)

Female

-0.04 (-0.29, 0.2)

-0.06 (-0.27, 0.15)

Male

-0.23 (-0.46, -0.01)

-0.31 (-0.55, -0.07)

BDE-47

All

0.11 (-0.06, 0.28)

0.16 (-0.01, 0.32)

Female

0.02 (-0.2, 0.25)

0 (-0.19, 0.18)

Male

0.22 (-0.04, 0.48)

0.31 (0.04, 0.58)

BDE-99

All

0.09 (-0.07, 0.25)

0.15 (-0.01, 0.3)

Female

0.1 (-0.1, 0.31)

0.13 (-0.04, 0.3)

Male

0.09 (-0.16, 0.34)

0.14 (-0.12, 0.41)

BDE-100

All

0.07 (-0.1, 0.23)

0.15 (-0.01, 0.31)

Female

0.03 (-0.2, 0.26)

0.05 (-0.14, 0.25)

Male

0.08 (-0.17, 0.34)

0.17 (-0.1, 0.44)

Table S3

Linear regression associations between natural log transformed polybrominated diphenyl ethers (PBDEs), measured in serum samples obtained during the first trimester, and internalizing and externalizing behaviors averaged across ages 1-5 years among Black mother-child pairs in metropolitan Atlanta, Georgia. Sensitivity analysis additionally adjusted for total lipids.

# Assign predictors, outcomes
predictors <- c("PBDE47.log", "PBDE99.log", "PBDE100.log")
outcomes <- c("Internalizing", "Externalizing")

# Loop regression
models <- list()
for (i in 1:length(outcomes)) {
  for (j in 1:length(predictors)) {
    model_name <- paste0(outcomes[i], " ~ ", predictors[j])
    formula <- as.formula(paste(outcomes[i], " ~ ", predictors[j], " + `Mean Child Age`  + `Infant Sex` + `Marital Status` + `Pre-pregnancy BMI` + `Maternal Age` + `Maternal Education` + Parity + `Substance Use` + total_lipids"))
    model <- lm(formula = formula, data = mix)
    models[[model_name]] <- tidy(model, conf.int = TRUE) %>%
      mutate(N = nobs(model)) %>%
      mutate(outcome = outcomes[i])
  }
}
chemicals_adj_lipid <- do.call(rbind, models) 

# Format 
reg_lipid <- chemicals_adj_lipid %>% 
  filter(term %in% c("PFHXS.log", "PFOS.log", "PFOA.log", "PFNA.log", "PBDE47.log", "PBDE99.log", "PBDE100.log"))  %>%
  dplyr::select(outcome, N, term, estimate, conf.low, conf.high)

# Combine regression results
# Table
table_s3 <- reg_lipid %>%
  mutate_if(is.numeric, round, digits = 2) %>%
  mutate(`Beta CI` = paste0(estimate, " (", conf.low, ", ", conf.high, ")"))  %>%
 mutate(term = recode(term, 
                     "PFHXS.log" = "PFHXS",
                     "PFOA.log" = "PFOA",
                     "PFNA.log" = "PFNA",
                     "PFOS.log" = "PFOS",
                     "PBDE47.log" = "BDE-47",
                     "PBDE99.log" = "BDE-99",
                     "PBDE100.log" = "BDE-100")) %>%
  dplyr::select(term, N, outcome, `Beta CI`) %>%
  pivot_wider(names_from = outcome, values_from = c("N", "Beta CI"),names_glue = "{outcome}_{.value}") %>%
  dplyr::select(term, Internalizing_N, `Internalizing_Beta CI`, Externalizing_N, `Externalizing_Beta CI`) %>%
  flextable() %>%
  separate_header() %>%
  font(fontname = "Tahoma") %>%
  set_header_labels(term = "Chemical") %>%
  autofit() %>%
  theme_alafoli()

table_s3

term

Internalizing

Externalizing

N

Beta CI

N

Beta CI

BDE-47

154

0.12 (-0.05, 0.29)

154

0.16 (-0.01, 0.33)

BDE-99

154

0.09 (-0.07, 0.25)

154

0.14 (-0.02, 0.3)

BDE-100

154

0.09 (-0.08, 0.25)

154

0.17 (0, 0.33)

Table S4

Linear regression associations between natural log transformed per- and polyfluoroalkyl substances (PFAS) and polybrominated diphenyl ethers (PBDEs), measured in serum samples obtained during the first trimester, and internalizing and externalizing behaviors averaged across ages 1-5 years among Black mother-child pairs in metropolitan Atlanta, Georgia. Sensitivity analysis removing preterm births.

# Removing preterm brith 
mix_ptb <- mix %>% filter(ptb == 0)
# Assign predictors, outcomes
predictors <- c("PFHXS.log", "PFOS.log", "PFOA.log", "PFNA.log", "PBDE47.log", "PBDE99.log", "PBDE100.log")
outcomes <- c("Internalizing", "Externalizing")

# Loop regression
models <- list()
for (i in 1:length(outcomes)) {
  for (j in 1:length(predictors)) {
    model_name <- paste0(outcomes[i], " ~ ", predictors[j])
    formula <- as.formula(paste(outcomes[i], " ~ ", predictors[j], " + `Mean Child Age`  + `Infant Sex` + `Marital Status` + `Pre-pregnancy BMI` + `Maternal Age` + `Maternal Education` + Parity + `Substance Use`"))
    model <- lm(formula = formula, data = mix_ptb)
    models[[model_name]] <- tidy(model, conf.int = TRUE) %>%
      mutate(outcome = outcomes[i])
  }
}
chemicals_adj_ptb <- do.call(rbind, models) 

# Format 
reg_ptb <- chemicals_adj_ptb %>% 
  filter(term %in% c("PFHXS.log", "PFOS.log", "PFOA.log", "PFNA.log", "PBDE47.log", "PBDE99.log", "PBDE100.log"))  %>%
  dplyr::select(outcome, term, estimate, conf.low, conf.high)

# Combine regression results
# Table
table_s4 <- reg_ptb %>%
  mutate_if(is.numeric, round, digits = 2) %>%
  mutate(Beta_CI = paste0(estimate, " (", conf.low, ", ", conf.high, ")"))  %>%
 mutate(term = recode(term, 
                     "PFHXS.log" = "PFHXS",
                     "PFOA.log" = "PFOA",
                     "PFNA.log" = "PFNA",
                     "PFOS.log" = "PFOS",
                     "PBDE47.log" = "BDE-47",
                     "PBDE99.log" = "BDE-99",
                     "PBDE100.log" = "BDE-100")) %>%
  dplyr::select(term, outcome, Beta_CI) %>%
  pivot_wider(names_from = outcome, values_from = (Beta_CI)) %>%
  group_by(term) %>%
  summarise(
    Internalizing = Internalizing,
    Externalizing = Externalizing) %>%
  flextable() %>%
  add_body_row(values = c("","Beta (95% CI)", "Beta (95% CI)")) %>%
  font(fontname = "Tahoma") %>%
  set_header_labels(term = "Chemical") %>%
  autofit() %>%
  theme_alafoli()

table_s4

Chemical

Internalizing

Externalizing

Beta (95% CI)

Beta (95% CI)

BDE-100

0.05 (-0.13, 0.22)

0.13 (-0.04, 0.31)

BDE-47

0.17 (-0.02, 0.35)

0.16 (-0.03, 0.35)

BDE-99

0.13 (-0.04, 0.3)

0.15 (-0.03, 0.32)

PFHXS

-0.18 (-0.35, -0.01)

-0.24 (-0.41, -0.07)

PFNA

-0.16 (-0.33, 0.01)

-0.22 (-0.4, -0.05)

PFOA

-0.11 (-0.29, 0.06)

-0.2 (-0.37, -0.03)

PFOS

-0.12 (-0.29, 0.05)

-0.24 (-0.41, -0.07)

Mixture Methods

Quantile G-Computation

Total Sample

# Define exposure, outcome, dataset
mixture <- mix %>% dplyr::select(PFHXS.log, PFOS.log, PFOA.log, PFNA.log, PBDE47.log, PBDE99.log, PBDE100.log, Internalizing, Externalizing, `Mean Child Age`, `Infant Sex`, `Maternal Age`, `Maternal Education`, Parity, `Pre-pregnancy BMI`, `Substance Use`)

exposure.names<-c("PFHXS.log", "PFOS.log", "PFOA.log","PFNA.log", "PBDE47.log", "PBDE99.log", "PBDE100.log")

outcomes <- c("Internalizing", "Externalizing")

# Build function for  qgcomp 

qgcompdata <- function(outcome) {
  model <- qgcomp.glm.noboot(as.formula(paste(outcome, "~", "PFHXS.log + PFOS.log + PFOA.log + PFNA.log + PBDE47.log + PBDE99.log + PBDE100.log + `Mean Child Age` + `Infant Sex` + `Maternal Age` +`Maternal Education` + Parity +`Pre-pregnancy BMI` + +`Substance Use`")),
                           expnms = exposure.names,
                           dat = mixture,
                           family = gaussian(),
                           q = 4)

 coefs <- data.frame(
    estimate = round(model$psi,2),
    lower_ci = round(model$ci.coef[-1, 1],2),
    upper_ci = round(model$ci.coef[-1, 2],2),
    outcome = outcome, 
    row.names = NULL
  )
adj_df <- data.frame()
adj_df <- rbind(adj_df, coefs)
print(adj_df)
}

# Run
results_all_int <- qgcompdata(outcome = "Internalizing")
results_all_ext <- qgcompdata(outcome = "Externalizing")

# Combine, format
qgcomp_all_df <- rbind(results_all_int, results_all_ext)
qgcomp_all_df <- qgcomp_all_df %>% mutate(group = "Total Mixture")

Chemical Specific

##### PFAS

# Define exposure, outcome dataset
mixture <- mix %>% dplyr::select(PFHXS.log, PFOS.log, PFOA.log, PFNA.log, PBDE47.log, PBDE99.log, PBDE100.log, Internalizing, Externalizing, `Mean Child Age`, `Infant Sex`, `Maternal Age`, `Maternal Education`, Parity, `Pre-pregnancy BMI`, `Insurance Type`, `Substance Use`, `Marital Status`)

outcomes <- c("Internalizing", "Externalizing")

exposure.names<-c("PFHXS.log", "PFOS.log", "PFOA.log", "PFNA.log")

# Function for  qgcomp 

qgcompdata <- function(outcome) {
  model <- qgcomp.glm.noboot(as.formula(paste(outcome, "~", "PFHXS.log + PFOS.log + PFOA.log + PFNA.log + `Mean Child Age` + `Infant Sex` + `Maternal Age` +`Maternal Education` + Parity +`Pre-pregnancy BMI` + `Insurance Type` +`Substance Use` +`Marital Status`")),
                  expnms = exposure.names,
                  dat = mixture,
                 family = gaussian(),
                 q = 4)

 coefs <- data.frame(
    estimate = round(model$psi,2),
    lower_ci = round(model$ci.coef[-1, 1],2),
    upper_ci = round(model$ci.coef[-1, 2],2),
    outcome = outcome, 
    row.names = NULL
  )
adj_df <- data.frame()
PFAS_df <- rbind(adj_df, coefs)
print(PFAS_df)
}

# Run
results_pfas_int <- qgcompdata(outcome = "Internalizing")
results_pfas_ext <- qgcompdata(outcome = "Externalizing")

qgcomp_pfas_df <- rbind(results_pfas_int, results_pfas_ext)

#### PBDE

# Define exposure, outcome, dataset

mixture <- mix %>% dplyr::select(PFHXS.log, PFOS.log, PFOA.log, PFNA.log, PBDE47.log, PBDE99.log, PBDE100.log, Internalizing, Externalizing, `Mean Child Age`, `Infant Sex`, `Maternal Age`, `Maternal Education`, Parity, `Pre-pregnancy BMI`, `Insurance Type`, `Substance Use`, `Marital Status`)

exposure.names<-c("PBDE47.log", "PBDE99.log", "PBDE100.log")

outcomes <- c("Internalizing", "Externalizing")

qgcompdata <- function(outcome) {
  model <- qgcomp.glm.noboot(as.formula(paste(outcome, "~", "PBDE47.log + PBDE99.log + PBDE100.log + `Mean Child Age` + `Infant Sex` + `Maternal Age` +`Maternal Education` + Parity +`Pre-pregnancy BMI` + `Insurance Type` + +`Substance Use` +`Marital Status`")),
                  expnms = exposure.names,
                  dat = mixture,
                 family = gaussian(),
                 q = 4)

 coefs <- data.frame(
    estimate = round(model$psi,2),
    lower_ci = round(model$ci.coef[-1, 1],2),
    upper_ci = round(model$ci.coef[-1, 2],2),
    outcome = outcome, 
    row.names = NULL
  )
adj_df <- data.frame()
adj_df <- rbind(adj_df, coefs)
print(adj_df)
}

# Run
results_pbde_int <- qgcompdata(outcome = "Internalizing")
results_pbde_ext <- qgcompdata(outcome = "Externalizing")

# Combine results
qgcomp_pbde_df <- rbind(results_pbde_int, results_pbde_ext)
qgcomp_chemicals_df <- rbind(qgcomp_pbde_df, qgcomp_pfas_df)
qgcomp_chemicals_df <-  qgcomp_chemicals_df %>% mutate(group = c("PBDE", "PBDE", "PFAS", "PFAS"))

Figure 2

Quantile g-computation cumulative effect (95% confidence intervals) and weights reflecting the associations between a mixture of natural log transformed per- and polyfluoroalkyl substances (PFAS, ng/mL) and polybrominated diphenyl ethers (PBDEs, pg/mL), measured in serum samples obtained during the first trimester, and internalizing and externalizing behaviors averaged across ages 1-5 years among Black mother-child pairs in metropolitan Atlanta, Georgia.

Models are adjusted for infant sex, mean child age, maternal age, maternal education, marital status, early pregnancy BMI, parity, and substance use. multiply by negative 1

##### TOTAL SAMPLE
# Forest plot
qgcomp_all_df$outcome <- factor(qgcomp_all_df$outcome, levels = c("Internalizing", "Externalizing"))
forest_all_qgcomp <- ggplot(qgcomp_all_df, aes(estimate, outcome)) +
  geom_pointrange(aes(xmin = lower_ci, xmax = upper_ci),
                  size = 2, linewidth = 2) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "grey25") +
  coord_flip() +
  scale_x_continuous(limits = c(-0.6, 0.4), breaks = c(-0.6, -0.4, -0.2, 0, 0.2, 0.4)) +
  ggtitle("Total Sample") + 
  theme_classic(base_family = "AppleGothic") +
  labs(x = "Beta (95 % Confidence Interval)", y = "CBCL Outcome") +
  theme(text = element_text(size = 25, color = "grey25", family = "Helvetica"),
        legend.title = element_text(color = "grey25"),
        axis.text = element_text(color = "grey25", size = 25),
        plot.title = element_text(hjust = 0.5))

# Weights
# Internalizing
exposure.names<-c("PFHXS.log", "PFOS.log", "PFOA.log","PFNA.log", "PBDE47.log", "PBDE99.log", "PBDE100.log")

model <- qgcomp.glm.noboot(Internalizing ~ PFHXS.log + PFOS.log + PFOA.log + PFNA.log + PBDE47.log + PBDE99.log + PBDE100.log + `Mean Child Age` + `Infant Sex` + `Maternal Age` +`Maternal Education` + Parity +`Pre-pregnancy BMI` +`Substance Use`,
                           expnms = exposure.names,
                           dat = mixture,
                           family = gaussian(),
                           q = 4)
# Extract weights
pos <- as.data.frame(model$pos.weights) 
# Format
pos <- pos %>% rename(weight = `model$pos.weights`) %>%
  mutate(chemical = c("BDE-47", "PFOS", "BDE-99", "PFOA"), outcome = "Internalizing", Direction = "Positive")
neg <- as.data.frame(model$neg.weights)
neg <- neg %>% rename(weight = `model$neg.weights`) %>%
  mutate(chemical = c("PFHXS", "PFNA", "BDE-100"), outcome = "Internalizing", Direction = "Negative")

neg$weight <- neg$weight*-1


int_all_weights_qgcomp <- rbind(pos, neg)



# Externalizing
model <- qgcomp.glm.noboot(Externalizing ~ PFHXS.log + PFOS.log + PFOA.log + PFNA.log + PBDE47.log + PBDE99.log + PBDE100.log + `Mean Child Age` + `Infant Sex` + `Maternal Age` +`Maternal Education` + Parity +`Pre-pregnancy BMI` +`Substance Use`,
                           expnms = exposure.names,
                           dat = mixture,
                           family = gaussian(),
                           q = 4)

# Extract weights
pos <- as.data.frame(model$pos.weights) 
pos <- pos %>% rename(weight = `model$pos.weights`) %>%
  mutate(chemical = c("BDE-47", "PFOA", "BDE-100", "BDE-99", "PFOS"), outcome = "Externalizing", Direction = "Positive")
neg <- as.data.frame(model$neg.weights)
neg <- neg %>% rename(weight = `model$neg.weights`) %>%
  mutate(chemical = c("PFHXS", "PFNA"), outcome = "Externalizing", Direction = "Negative")

neg$weight <- neg$weight*-1

ext_all_weights_qgcomp <- rbind(pos, neg)

# Plot
int_all_weights_qgcomp <- ggplot(int_all_weights_qgcomp, aes(fct_reorder(chemical, weight), weight, fill = Direction)) +
  geom_col() +
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey25") +
  labs(x = "Chemical", y = "Weight", fill = "Direction", title = "Internalizing") +
scale_fill_manual(values = c("salmon3", "steelblue3")) +
  scale_y_continuous(limits = c(-1, 1), breaks = c(-1, -0.5, 0, 0.5, 1)) +
  theme_classic(base_family = "AppleGothic") +
  theme(text = element_text(size = 25, color = "grey25"),
        legend.title = element_text(color = "grey25"),
        axis.text.x = element_text(color = "grey25", angle = 45, hjust = 1, size = 25),
        axis.text.y = element_text(color = "grey25", size = 25),
        legend.position = "none",
        plot.title = element_text(hjust = 0.5))

ext_all_weights_qgcomp <- ggplot(ext_all_weights_qgcomp, aes(fct_reorder(chemical, weight), weight, fill = Direction)) +
  geom_col() +
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey25") +
  labs(x = "Chemical", y = "Weight", fill = "Direction", title = "Externalizing") +
    scale_y_continuous(limits = c(-1, 1), breaks = c(-1, -0.5, 0, 0.5, 1)) +
scale_fill_manual(values = c("salmon3", "steelblue3")) +
  theme_classic(base_family = "AppleGothic") +
  theme(text = element_text(size = 25, color = "grey25"),
        legend.title = element_text(color = "grey25"),
        axis.text.x = element_text(color = "grey25", angle = 45, hjust = 1, size = 25),
        axis.text.y = element_text(color = "grey25", size = 25),
        legend.position = "bottom",
        plot.title = element_text(hjust = 0.5),
        legend.key.size = unit(3, "cm"),
        legend.text = element_text(size = 25))

##### CHEMCIAL SPECIFIC

# Forest plot
qgcomp_pfas_df$outcome <- factor(qgcomp_pfas_df$outcome, levels = c("Internalizing", "Externalizing"))
forest_pfas_qgcomp <- ggplot(qgcomp_pfas_df, aes(estimate, outcome)) +
  geom_pointrange(aes(xmin = lower_ci, xmax = upper_ci),
                  size = 2, linewidth = 2, color = "black", position = position_dodge(width = 0.25)) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "grey25") +
  scale_x_continuous(limits = c(-0.6, 0.4), breaks = c(-0.6, -0.4, -0.2, 0, 0.2, 0.4)) +
  ggtitle("PFAS") + 
  coord_flip() +
  theme_classic(base_family = "AppleGothic") +
  labs(x = "Beta (95 % Confidence Interval)", y = "CBCL Outcome", linetype = "Chemical Class") +
  theme(text = element_text(size = 25, color = "grey25", family = "Helvetica"),
        legend.title = element_text(color = "grey25"),
        axis.text = element_text(color = "grey25", size = 25),
        plot.title = element_text(hjust = 0.5))

forest_pbde_qgcomp <- ggplot(qgcomp_pbde_df, aes(estimate, outcome)) +
  geom_pointrange(aes(xmin = lower_ci, xmax = upper_ci),
                  size = 2, linewidth = 2, color = "black", position = position_dodge(width = 0.25)) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "grey25") +
  scale_x_continuous(limits = c(-0.6, 0.4), breaks = c(-0.6, -0.4, -0.2, 0, 0.2, 0.4)) +
  ggtitle("PBDE") +
  coord_flip() +
  theme_classic(base_family = "AppleGothic") +
  labs(x = "Beta (95 % Confidence Interval)", y = "CBCL Outcome", linetype = "Chemical Class") +
  theme(text = element_text(size = 25, color = "grey25", family = "Helvetica"),
        legend.title = element_text(color = "grey25"),
        axis.text = element_text(color = "grey25", size = 25),
        plot.title = element_text(hjust = 0.5))

# Weights
# PFAS 
exposure.names<-c("PFHXS.log", "PFOS.log", "PFOA.log", "PFNA.log")

model <- qgcomp.glm.noboot(Internalizing ~ PFHXS.log + PFOS.log + PFOA.log + PFNA.log + `Mean Child Age` + `Infant Sex` + `Maternal Age` +`Maternal Education` + Parity +`Pre-pregnancy BMI` +`Substance Use`,
                           expnms = exposure.names,
                           dat = mixture,
                           family = gaussian(),
                           q = 4)
# Extract weights
pos <- as.data.frame(model$pos.weights) 
pos <- pos %>% rename(weight = `model$pos.weights`) %>%
  mutate(chemical = c("PFOS", "PFOA"), outcome = "Internalizing", Direction = "Positive")
neg <- as.data.frame(model$neg.weights)
neg <- neg %>% rename(weight = `model$neg.weights`) %>%
  mutate(chemical = c("PFNA", "PFHXS"), outcome = "Internalizing", Direction = "Negative")

neg$weight <- neg$weight*-1

weights_int_pfas_qgcomp <- rbind(pos, neg)

# Externalizing
model <- qgcomp.glm.noboot(Externalizing ~ PFHXS.log + PFOS.log + PFOA.log + PFNA.log + `Mean Child Age` + `Infant Sex` + `Maternal Age` +`Maternal Education` + Parity +`Pre-pregnancy BMI` +`Substance Use`,
                           expnms = exposure.names,
                           dat = mixture,
                           family = gaussian(),
                           q = 4)
# Extract weights
pos <- as.data.frame(model$pos.weights) 
pos <- pos %>% rename(weight = `model$pos.weights`) %>%
  mutate(chemical = c("PFOA", "PFOS"), outcome = "Externalizing", Direction = "Positive")
neg <- as.data.frame(model$neg.weights)
neg <- neg %>% rename(weight = `model$neg.weights`) %>%
  mutate(chemical = c("PFNA", "PFHXS"), outcome = "Externalizing", Direction = "Negative")

neg$weight <- neg$weight*-1

weights_ext_pfas_qgcomp <- rbind(pos, neg)

# PBDE
exposure.names<-c("PBDE47.log", "PBDE99.log", "PBDE100.log")

model <- qgcomp.glm.noboot(Internalizing ~ PBDE47.log + PBDE99.log + PBDE100.log + `Mean Child Age` + `Infant Sex` + `Maternal Age` +`Maternal Education` + Parity +`Pre-pregnancy BMI` +`Substance Use`,
                           expnms = exposure.names,
                           dat = mixture,
                           family = gaussian(),
                           q = 4)
pos <- as.data.frame(model$pos.weights) 
pos <- pos %>% rename(weight = `model$pos.weights`) %>%
  mutate(chemical = c("BDE-47", "BDE-99"), outcome = "Internalizing", Direction = "Positive")
neg <- as.data.frame(model$neg.weights)
neg <- neg %>% rename(weight = `model$neg.weights`) %>%
  mutate(chemical = c("BDE-100"), outcome = "Internalizing", Direction = "Negative")

neg$weight <- neg$weight*-1

# Extract weights
weights_int_pbde_qgcomp <- rbind(pos, neg)

# Externalizing
model <- qgcomp.glm.noboot(Externalizing ~ PBDE47.log + PBDE99.log + PBDE100.log + `Mean Child Age` + `Infant Sex` + `Maternal Age` +`Maternal Education` + Parity +`Pre-pregnancy BMI` +`Substance Use`,
                           expnms = exposure.names,
                           dat = mixture,
                           family = gaussian(),
                           q = 4)
# Extract weights
pos <- as.data.frame(model$pos.weights) 
pos <- pos %>% rename(weight = `model$pos.weights`) %>%
  mutate(chemical = c("BDE-47", "BDE-99", "BDE-100"), outcome = "Externalizing", Direction = "Positive")

weights_ext_pbde_qgcomp <- pos

# Plot
int_pfas_weights_qgcomp <- ggplot(weights_int_pfas_qgcomp, aes(fct_reorder(chemical, weight), weight, fill = Direction)) +
  geom_col() +
  geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
  labs(x = "Chemical", y = "Weight", fill = "Direction", title = "Internalizing") +
scale_fill_manual(values = c("salmon3", "steelblue3")) +
  scale_y_continuous(limits = c(-1, 1), breaks = c(-1, -0.5, 0, 0.5, 1)) +
  theme_classic(base_family = "AppleGothic") +
  theme(text = element_text(size = 20, color = "grey25"),
        legend.title = element_text(color = "grey25"),
        axis.text.x = element_text(color = "grey25", angle = 45, hjust = 1, size = 25),
        axis.text.y = element_text(color = "grey25", size = 25),
        legend.position = "none",
        plot.title = element_text(hjust = 0.5))

ext_pfas_weights_qgcomp <- ggplot(weights_ext_pfas_qgcomp, aes(fct_reorder(chemical, weight), weight, fill = Direction)) +
  geom_col() +
  geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
  labs(x = "Chemical", y = "Weight", fill = "Direction", title = "Externalizing") +
scale_fill_manual(values = c("salmon3", "steelblue3")) +
    scale_y_continuous(limits = c(-1, 1), breaks = c(-1, -0.5, 0, 0.5, 1)) +
  theme_classic(base_family = "AppleGothic") +
  theme(text = element_text(size = 25, color = "grey25"),
        legend.title = element_text(color = "grey25"),
        axis.text.x = element_text(color = "grey25", angle = 45, hjust = 1, size = 25),
        axis.text.y = element_text(color = "grey25", size = 25),
        legend.position = "bottom",
        plot.title = element_text(hjust = 0.5),
        legend.key.size = unit(3, "cm"),
        legend.text = element_text(size = 25))


int_pbde_weights_qgcomp <- ggplot(weights_int_pbde_qgcomp, aes(fct_reorder(chemical, weight), weight, fill = Direction)) +
  geom_col() +
  geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
  labs(x = "Chemical", y = "Weight", fill = "Direction", title = "Internalizing") +
scale_fill_manual(values = c("salmon3", "steelblue3")) +
  theme_classic(base_family = "AppleGothic") +
  scale_y_continuous(limits = c(-1, 1), breaks = c(-1, -0.5, 0, 0.5, 1)) +
  theme(text = element_text(size = 25, color = "grey25"),
        legend.title = element_text(color = "grey25"),
        axis.text.x = element_text(color = "grey25", angle = 45, hjust = 1, size = 25),
        axis.text.y = element_text(color = "grey25", size = 25),
        legend.position = "none",
        plot.title = element_text(hjust = 0.5))

ext_pbde_weights_qgcomp <- ggplot(weights_ext_pbde_qgcomp, aes(fct_reorder(chemical, weight), weight, fill = Direction)) +
  geom_col() +
  geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
  labs(x = "Chemical", y = "Weight", fill = "Direction", title = "Externalizing") +
scale_fill_manual(values = c("steelblue3")) +
  theme_classic(base_family = "AppleGothic") +
  scale_y_continuous(limits = c(-1, 1), breaks = c(-1, -0.5, 0, 0.5, 1)) +
  theme(text = element_text(size = 25, color = "grey25"),
        legend.title = element_text(color = "grey25"),
        axis.text.x = element_text(color = "grey25", angle = 45, hjust = 1, size = 25),
        axis.text.y = element_text(color = "grey25", size = 25),
        legend.position = "bottom",
        plot.title = element_text(hjust = 0.5),
        legend.key.size = unit(3, "cm"),
        legend.text = element_text(size = 25))

figure_3 <- plot_grid(forest_all_qgcomp, forest_pfas_qgcomp, forest_pbde_qgcomp, int_all_weights_qgcomp,
int_pfas_weights_qgcomp, int_pbde_weights_qgcomp,ext_all_weights_qgcomp, ext_pfas_weights_qgcomp, 
ext_pbde_weights_qgcomp, ncol = 3, nrow = 3, align = "v")

figure_3

Table S5

Association between the total mixture of per- and polyfluoroalkyl substances (PFAS, ng/mL) and polybrominated diphenyl ethers (PBDEs, pg/mL), measured in first trimester serum samples, a mixture of PFAS alone, and a mixture of PBDEs alone, estimated using quantile g-computation, and internalizing and externalizing behaviors averaged across ages 1-5 years among Black mother-child pairs in metropolitan Atlanta, Georgia.

Models are adjusted for infant sex, mean child age, maternal age, maternal education, marital status, early pregnancy BMI, parity, and substance use.

# Combine results
table_s6 <- rbind(qgcomp_all_df, qgcomp_chemicals_df)

# Table
table_s6 <- table_s6 %>%
  mutate_if(is.numeric, round, digits = 2) %>%
  mutate("Beta (CI)" = paste0(estimate, " (", lower_ci, ", ", upper_ci, ")"))  %>%
  dplyr::select(group, outcome, `Beta (CI)`) %>%
  pivot_wider(names_from = outcome, values_from = `Beta (CI)`) %>%
  flextable() %>%
  font(fontname = "Tahoma") %>%
  add_body_row(values = c("","Beta (95% CI)", "Beta (95% CI)")) %>%
  merge_v(j = "group") %>%
  set_header_labels(
  group = "") %>%
  autofit() %>%
  theme_alafoli()

table_s6

Internalizing

Externalizing

Beta (95% CI)

Beta (95% CI)

Total Mixture

-0.04 (-0.28, 0.2)

-0.06 (-0.3, 0.18)

PBDE

0.13 (-0.03, 0.3)

0.2 (0.04, 0.36)

PFAS

-0.2 (-0.39, -0.01)

-0.28 (-0.46, -0.1)

Bayeasian Kernal Machine Regression

Total Sample

# Define dataset
mixture <- mix %>% dplyr::select(PFHXS.log, PFOS.log, PFOA.log, PFNA.log, PBDE47.log, PBDE99.log, PBDE100.log, Internalizing, Externalizing, `Mean Child Age`, `Infant Sex`, `Maternal Age`, `Maternal Education`, Parity, ppbmi, `Substance Use`)

# Define mixture, covariates
df.exp <- mixture %>% dplyr::select(PFHXS.log, PFOS.log, PFOA.log, PFNA.log, PBDE47.log, PBDE99.log, PBDE100.log)

df.cov <- mixture %>% dplyr::select(`Mean Child Age`, `Infant Sex`, `Maternal Age`, `Maternal Education`, Parity, ppbmi, `Substance Use`)

df.cov <- sapply(df.cov, as.numeric)

# BKMR


# Internalizing
set.seed(123)
fit_int <- kmbayes(y = mixture$Internalizing,
                Z = df.exp,
                X = df.cov,
                iter = 10000, 
                verbose = FALSE, 
                varsel = TRUE,
                family='gaussian')

# Extract pips
pips_int <- ExtractPIPs(fit_int)
pips_int <- pips_int %>% mutate(outcome = "Internalizing", group = "Total Sample")

# Externalizing
set.seed(123)
fit_ext <- kmbayes(y = mixture$Externalizing,
                Z = df.exp,
                X = df.cov,
                iter = 10000, 
                verbose = FALSE, 
                varsel = TRUE,
                family='gaussian')

# Extract pips
pips_ext <- ExtractPIPs(fit_ext)
pips_ext <- pips_ext %>% mutate(outcome = "Externalizing", group = "Total Sample")

Chemical Specific

## PFAS 
# Define dataset
mixture <- mix %>% dplyr::select(PFHXS.log, PFOS.log, PFOA.log, PFNA.log, PBDE47.log, PBDE99.log, PBDE100.log, Internalizing, Externalizing, `Mean Child Age`, `Infant Sex`, `Maternal Age`, `Maternal Education`, Parity, ppbmi, `Substance Use`)

# Define mixture, covariates
df.exp <- mixture %>% dplyr::select(PFHXS.log, PFOS.log, PFOA.log, PFNA.log)

df.cov <- mixture %>% dplyr::select(`Mean Child Age`, `Infant Sex`, `Maternal Age`, `Maternal Education`, Parity, ppbmi,`Substance Use`)

df.cov <- sapply(df.cov, as.numeric)

# BKMR
# Internalizing
set.seed(123)
fit_int_pfas <- kmbayes(y = mixture$Internalizing,
                Z = df.exp,
                X = df.cov,
                iter = 10000, 
                verbose = FALSE, 
                varsel = TRUE,
                family='gaussian')


# Extract pips
pips_int_pfas <- ExtractPIPs(fit_int_pfas)
pips_int_pfas <- pips_int_pfas %>% mutate(outcome = "Internalizing", group = "PFAS")


# Externalizing
set.seed(123)
fit_ext_pfas <- kmbayes(y = mixture$Externalizing,
                Z = df.exp,
                X = df.cov,
                iter = 10000, 
                verbose = FALSE, 
                varsel = TRUE,
                family='gaussian')

# Extract pips
pips_ext_pfas <- ExtractPIPs(fit_ext_pfas)
pips_ext_pfas <- pips_ext_pfas %>% mutate(outcome = "Externalizing", group = "PFAS")


## PBDE
# Define dataset
mixture <- mix %>% dplyr::select(PFHXS.log, PFOS.log, PFOA.log, PFNA.log, PBDE47.log, PBDE99.log, PBDE100.log, Internalizing, Externalizing, `Mean Child Age`, `Infant Sex`, `Maternal Age`, `Maternal Education`, Parity, ppbmi, `Substance Use`) 

# Define mixture, covariates
df.exp <- mixture %>% dplyr::select(PBDE47.log, PBDE99.log, PBDE100.log)

df.cov <- mixture %>% dplyr::select(`Mean Child Age`, `Infant Sex`, `Maternal Age`, `Maternal Education`, Parity, ppbmi, `Substance Use`)

df.cov <- sapply(df.cov, as.numeric)

# BKMR


# Internalizing 
set.seed(123)
fit_int_pbde <- kmbayes(y = mixture$Internalizing,
                Z = df.exp,
                X = df.cov,
                iter = 10000, 
                verbose = FALSE, 
                varsel = TRUE,
                family='gaussian')

# Extract pips
pips_int_pbde <- ExtractPIPs(fit_int_pbde)
pips_int_pbde <- pips_int_pbde %>% mutate(outcome = "Internalizing", group = "PBDE")

# Externalizing
set.seed(123)
fit_ext_pbde <- kmbayes(y = mixture$Externalizing,
                Z = df.exp,
                X = df.cov,
                iter = 10000, 
                verbose = FALSE, 
                varsel = TRUE,
                family='gaussian')

# Extract pips
pips_ext_pbde <- ExtractPIPs(fit_ext_pbde)
pips_ext_pbde <- pips_ext_pbde %>% mutate(outcome = "Externalizing", group = "PBDE")

Table S6

Posterior inclusion probabilities (PIPs) reflecting the importance of each exposure in Bayesian Kernal Machine Regression models estimating the association between natural log transformed per- and polyfluoroalkyl substances (PFAS, ng/mL) and polybrominated diphenyl ethers (PBDEs, pg/mL), measured in first trimester serum samples, a mixture of PFAS alone, and a mixture of PBDEs alone.

Models are adjusted for infant sex, mean child age, maternal age, maternal education, marital status, early pregnancy BMI, parity, and substance use.

# Combine results
bkmr <- rbind(pips_int, pips_ext, pips_int_pfas, pips_ext_pfas, pips_int_pbde, pips_ext_pbde)

# Table
table_s7 <- bkmr %>%
  mutate_if(is.numeric, round, digits = 2) %>%
  dplyr::select(outcome, group, variable, PIP) %>%
  mutate(variable = recode(variable, 
                     "PFHXS.log" = "PFHXS",
                     "PFOA.log" = "PFOA",
                     "PFNA.log" = "PFNA",
                     "PFOS.log" = "PFOS",
                     "PBDE47.log" = "BDE-47",
                     "PBDE99.log" = "BDE-99",
                     "PBDE100.log" = "BDE-100")) %>%
  pivot_wider(names_from = outcome, values_from = PIP) %>%
  as_grouped_data(groups = "group") %>%
  flextable() %>%
  add_body_row(values = c("", "", "PIP", "PIP")) %>%
  set_header_labels(group = "", variable = "Chemical") %>%
  font(fontname = "Tahoma") %>%
  autofit() %>%
  theme_alafoli()

table_s7

Chemical

Internalizing

Externalizing

PIP

PIP

Total Sample

PFHXS

0.47

0.44

PFOS

0.27

0.31

PFOA

0.24

0.16

PFNA

0.34

0.32

BDE-47

0.46

0.28

BDE-99

0.37

0.27

BDE-100

0.26

0.36

PFAS

PFHXS

0.50

0.48

PFOS

0.22

0.32

PFOA

0.19

0.23

PFNA

0.35

0.41

PBDE

BDE-47

0.54

0.44

BDE-99

0.56

0.46

BDE-100

0.30

0.38

Figure 3

Cumulative effect (95% credible intervals) of the BKMR mixture of natural log transformed per- and polyfluoroalkyl substances (PFAS) and polybrominated diphenyl ethers (PBDEs, pg/mL), measured in serum samples obtained during the first trimester, and internalizing and externalizing behaviors averaged across ages 1-5 years among Black mother-child pairs in metropolitan Atlanta, Georgia.

Models are adjusted for infant sex, mean child age, maternal age, maternal education, marital status, early pregnancy BMI, parity, and substance use.

### TOTAL SAMPLE
df.exp <- mixture %>% dplyr::select(PFHXS.log, PFOS.log, PFOA.log, PFNA.log, PBDE47.log, PBDE99.log, PBDE100.log)

# Internalizing 
risk.overall_int <- OverallRiskSummaries(fit =fit_int, y = mixture$Internalizing, Z = as.matrix(df.exp), X = as.matrix(df.cov), 
                                      qs = seq(0.25, 0.75, by = 0.05), 
                                      q.fixed = 0.5, method = "exact")

# Externalizing

risk.overall_ext <- OverallRiskSummaries(fit = fit_ext, y = mixture$Externalizing, Z = as.matrix(df.exp), X = as.matrix(df.cov), 
                                      qs = seq(0.25, 0.75, by = 0.05), 
                                      q.fixed = 0.5, method = "exact")

# Merge
risk.overall_int <- risk.overall_int %>% mutate(outcome = "Internalizing") 
risk.overall_ext <- risk.overall_ext  %>% mutate(outcome = "Externalizing")
risk.overall <- rbind(risk.overall_int, risk.overall_ext)


## PFAS
df.exp <- mixture %>% dplyr::select(PFHXS.log, PFOS.log, PFOA.log, PFNA.log)

# Internalizing
risk.overall_int_pfas <- OverallRiskSummaries(fit = fit_int_pfas, y = mixture$Internalizing, Z = as.matrix(df.exp), X = as.matrix(df.cov), 
                                      qs = seq(0.25, 0.75, by = 0.05), 
                                      q.fixed = 0.5, method = "exact")

# Externalizing

risk.overall_ext_pfas <- OverallRiskSummaries(fit = fit_ext_pfas, y = mixture$Externalizing, Z = as.matrix(df.exp), X = as.matrix(df.cov), 
                                      qs = seq(0.25, 0.75, by = 0.05), 
                                      q.fixed = 0.5, method = "exact")

# Merge
risk.overall_int_pfas <- risk.overall_int_pfas %>% mutate(outcome = "Internalizing") 
risk.overall_ext_pfas <- risk.overall_ext_pfas  %>% mutate(outcome = "Externalizing")
risk.overall_pfas <- rbind(risk.overall_int_pfas, risk.overall_ext_pfas)

## PBDE
df.exp <- mixture %>% dplyr::select(PBDE47.log, PBDE99.log, PBDE100.log)

# Internalizing
risk.overall_int_pbde <- OverallRiskSummaries(fit = fit_int_pbde, y = mixture$Internalizing, Z = as.matrix(df.exp), X = as.matrix(df.cov), 
                                      qs = seq(0.25, 0.75, by = 0.05), 
                                      q.fixed = 0.5, method = "exact")

# Externalizing

risk.overall_ext_pbde <- OverallRiskSummaries(fit = fit_ext_pbde, y = mixture$Externalizing, Z = as.matrix(df.exp), X = as.matrix(df.cov), 
                                      qs = seq(0.25, 0.75, by = 0.05), 
                                      q.fixed = 0.5, method = "exact")

# Merge
risk.overall_int_pbde <- risk.overall_int_pbde %>% mutate(outcome = "Internalizing") 
risk.overall_ext_pbde <- risk.overall_ext_pbde  %>% mutate(outcome = "Externalizing")
risk.overall_pbde <- rbind(risk.overall_int_pbde, risk.overall_ext_pbde)

# Plot

risk.overall_int_bkmr <- ggplot(risk.overall_int, aes(quantile, est, ymin = est - 1.96*sd, ymax = est + 1.96*sd)) + 
  geom_pointrange(size = 1.5, linewidth = 1.5, color = "grey25") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey25") +
  scale_y_continuous(limits = c(-0.4, 0.4), breaks = c(-0.4, -0.2, 0, 0.2, 0.4)) + 
  ggtitle("Total Sample
Internalizing") +
  theme_classic(base_family = "AppleGothic") +
  labs(x = "Quantile", y = "Beta (95 % Credible Interval)") +
  theme(text = element_text(size = 25, color = "grey25"),
        axis.text = element_text(color = "grey25", size = 25), 
        legend.title = element_blank(),
        plot.title = element_text(hjust = 0.5))
risk.overall_ext_bkmr <- ggplot(risk.overall_ext, aes(quantile, est, ymin = est - 1.96*sd, ymax = est + 1.96*sd)) + 
  geom_pointrange(size = 1.5, linewidth = 1.5, color =  "grey55") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey25") +
  scale_x_continuous(limits = c(0.25, 0.75), breaks = c(0.30, 0.40, 0.50, 0.60, 0.70)) + 
  scale_y_continuous(limits = c(-0.4, 0.4), breaks = c(-0.4, -0.2, 0, 0.2, 0.4)) + 
  theme_classic() +
  labs(x = "Quantile", y = "Beta (95 % Credible Interval)", title = "Externalizing") +
  theme(text = element_text(size = 25, color = "grey25"),
        axis.text = element_text(color = "grey25", size = 25), 
        legend.title = element_blank(),
         plot.title = element_text(hjust = 0.5))
risk.overall_int_pfas_bkmr <- ggplot(risk.overall_int_pfas, aes(quantile, est, ymin = est - 1.96*sd, ymax = est + 1.96*sd)) + 
  geom_pointrange(size = 1.5, linewidth = 1.5, color = "grey25") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey25") +
  scale_y_continuous(limits = c(-0.4, 0.4), breaks = c(-0.4, -0.2, 0, 0.2, 0.4)) + 
  ggtitle("PFAS
     Internalizing") +
  theme_classic(base_family = "AppleGothic") +
  labs(x = "Quantile", y = "Beta (95 % Credible Interval)") +
  theme(text = element_text(size = 25, color = "grey25"),
        axis.text = element_text(color = "grey25", size = 25), 
        legend.title = element_blank(),
        plot.title = element_text(hjust = 0.5))
risk.overall_ext_pfas_bkmr <- ggplot(risk.overall_ext_pfas, aes(quantile, est, ymin = est - 1.96*sd, ymax = est + 1.96*sd)) + 
  geom_pointrange(size = 1.5, linewidth = 1.5, color = "grey55") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey25") +
  scale_y_continuous(limits = c(-0.4, 0.4), breaks = c(-0.4, -0.2, 0, 0.2, 0.4)) + 
  theme_classic(base_family = "AppleGothic") +
  labs(x = "Quantile", y = "Beta (95 % Credible Interval)", title = "Externalizing") +
  theme(text = element_text(size = 25, color = "grey25"),
        axis.text = element_text(color = "grey25", size = 25), 
        legend.title = element_blank(),
        plot.title = element_text(hjust = 0.5))
risk.overall_int_pbde_bkmr <- ggplot(risk.overall_int_pbde, aes(quantile, est, ymin = est - 1.96*sd, ymax = est + 1.96*sd)) + 
  geom_pointrange(size = 1.5, linewidth = 1.5, color = "grey25") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey25") +
  scale_y_continuous(limits = c(-0.4, 0.4), breaks = c(-0.4, -0.2, 0, 0.2, 0.4)) + 
  ggtitle("PBDE
    Internalizing") +
  theme_classic(base_family = "AppleGothic") +
  labs(x = "Quantile", y = "Beta (95 % Credible Interval)") +
  theme(text = element_text(size = 25, color = "grey25"),
        axis.text = element_text(color = "grey25", size = 25), 
        legend.title = element_blank(),
        plot.title = element_text(hjust = 0.5))
risk.overall_ext_pbde_bkmr <- ggplot(risk.overall_ext_pbde, aes(quantile, est, ymin = est - 1.96*sd, ymax = est + 1.96*sd)) + 
  geom_pointrange(size = 1.5, linewidth = 1.5, color = "grey55") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey25") +
  scale_y_continuous(limits = c(-0.4, 0.4), breaks = c(-0.4, -0.2, 0, 0.2, 0.4)) + 
  theme_classic(base_family = "AppleGothic") +
  labs(x = "Quantile", y = "Beta (95 % Credible Interval)", title = "Externalizing") +
  theme(text = element_text(size = 25, color = "grey25"),
        axis.text = element_text(color = "grey25", size = 25), 
        legend.title = element_blank(),
        plot.title = element_text(hjust = 0.5))
figure_4 <- plot_grid(risk.overall_int_bkmr, risk.overall_int_pfas_bkmr, risk.overall_int_pbde_bkmr, risk.overall_ext_bkmr,risk.overall_ext_pfas_bkmr, risk.overall_ext_pbde_bkmr, ncol = 3, row = 2)
figure_4

Figure S2

Univariate exposure-response functions indicating change in internalizing and externalizing behaviors resulting from individual natural log transformed per- and polyfluoroalkyl substances (PFAS, ng/mL) and polybrominated diphenyl ethers (PBDEs, pg/mL), measured in first trimester serum samples, while fixing remaining exposures in the mixture at the 50th percentile, estimated using Bayesian Kernal Machine Regression.

Models are adjusted for infant sex, mean child age, maternal age, maternal education, marital status, early pregnancy BMI, parity, and substance use.

# Build labels
chemicals_label <- c("PFHXS", "PFOS", "PFOA", "PFNA", "BDE-47", "BDE-99", "BDE-100")
names(chemicals_label) <- c("PFHXS.log", "PFOS.log", "PFOA.log", "PFNA.log", "PBDE47.log", "PBDE99.log", "PBDE100.log")

# Internalizing
pred.resp.univar <- PredictorResponseUnivar(fit = fit_int)

int_univar <- ggplot(pred.resp.univar, aes(z, est, ymin = est - 1.96*se, ymax = est + 1.96*se)) + 
  geom_smooth(stat = "identity", color  ="grey25", size = 2) + 
  facet_wrap(~ variable, label = labeller(variable = chemicals_label)) +
  scale_y_continuous(breaks = c(-1, 0, 1)) + 
  ylab("Internalizing") +
  xlab("Chemical") +
  theme_bw(base_family = "AppleGothic") +
  theme(text = element_text(size = 30, color = "grey25"),
      axis.text = element_text(color = "grey25", size = 30),
        strip.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank())

# Externalizing
pred.resp.univar <- PredictorResponseUnivar(fit = fit_ext)


ext_univar <- ggplot(pred.resp.univar, aes(z, est, ymin = est - 1.96*se, ymax = est + 1.96*se)) + 
  geom_smooth(stat = "identity", color = "grey25", size = 2) + 
  facet_wrap(~ variable, label = labeller(variable = chemicals_label)) +
  scale_y_continuous(breaks = c(-1, 0, 1)) + 
  xlab("Chemical") +
  ylab("Externalizing") +
  theme_bw(base_family = "AppleGothic") +
  theme(text = element_text(size = 30, color = "grey25"),
      axis.text = element_text(color = "grey25", size = 30),
        strip.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank())

figure_s2 <- plot_grid(int_univar, ext_univar, ncol = 2, align = "h")
figure_s2

Figure S3

Bivariate exposure-response functions indicating change in internalizing and externalizing behaviors resulting from individual natural log transformed per- and polyfluoroalkyl substances (PFAS, ng/mL) and polybrominated diphenyl ethers (PBDEs, pg/mL), measured in first trimester serum samples, while fixing remaining exposures in the mixture at specific percentiles, estimated using Bayesian Kernal Machine Regression.

Models are adjusted for infant sex, mean child age, maternal age, maternal education, marital status, early pregnancy BMI, parity, and substance use.

# Define mixture
df.exp <- mixture %>% dplyr::select(PFHXS.log, PFOS.log, PFOA.log, PFNA.log, PBDE47.log, PBDE99.log, PBDE100.log)

# Internalizing

pred.resp.bivar <- PredictorResponseBivar(fit = fit_int, min.plot.dist = 1)
pred.resp.bivar.levels <- PredictorResponseBivarLevels(
  pred.resp.df = pred.resp.bivar, 
                                                       
  Z = as.matrix(df.exp), qs = c(0.25, 0.5, 0.75))

# Plot

# Create labels
chemicals_label <- c("PFHXS", "PFOS", "PFOA", "PFNA", "BDE-47", "BDE-99", "BDE-100")
names(chemicals_label) <- c("PFHXS.log", "PFOS.log", "PFOA.log", "PFNA.log", "PBDE47.log", "PBDE99.log", "PBDE100.log")

pred.resp.bivar.levels$variable1 <- factor(pred.resp.bivar.levels$variable1, levels = c("PFHXS.log", "PFOS.log", "PFOA.log", "PFNA.log", "PBDE47.log", "PBDE99.log", "PBDE100.log"))
pred.resp.bivar.levels$variable2 <- factor(pred.resp.bivar.levels$variable2, levels = c("PFHXS.log", "PFOS.log", "PFOA.log", "PFNA.log", "PBDE47.log", "PBDE99.log", "PBDE100.log"))
 

figure_s21 <- ggplot(pred.resp.bivar.levels, aes(z1, est)) + 
    geom_smooth(aes(col = quantile), stat = "identity") + 
    facet_grid(variable2 ~ variable1, label = labeller(variable1 = chemicals_label, variable2 = chemicals_label)) +
 scale_y_continuous(breaks = c(-0.4, 0, 0.4)) + 
  xlab("Chemical") +
  ylab("Internalizing") +
  scale_color_manual(values = c("steelblue3", "salmon3", "springgreen4")) +
  theme_bw(base_family = "AppleGothic") +
  theme(text = element_text(size = 30, color = "grey25"),
      axis.text = element_text(color = "grey25", size = 30),
        strip.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
       legend.position = "bottom")
# Externalizing
pred.resp.bivar <- PredictorResponseBivar(fit = fit_ext, min.plot.dist = 1)
pred.resp.bivar.levels <- PredictorResponseBivarLevels(
  pred.resp.df = pred.resp.bivar, 
                                                       
  Z = as.matrix(df.exp), qs = c(0.25, 0.5, 0.75))

pred.resp.bivar.levels$variable1 <- factor(pred.resp.bivar.levels$variable1, levels = c("PFHXS.log", "PFOS.log", "PFOA.log", "PFNA.log", "PBDE47.log", "PBDE99.log", "PBDE100.log"))
pred.resp.bivar.levels$variable2 <- factor(pred.resp.bivar.levels$variable2, levels = c("PFHXS.log", "PFOS.log", "PFOA.log", "PFNA.log", "PBDE47.log", "PBDE99.log", "PBDE100.log"))

# Plot

figure_s22 <- ggplot(pred.resp.bivar.levels, aes(z1, est)) + 
    geom_smooth(aes(col = quantile), stat = "identity") + 
    facet_grid(variable2 ~ variable1, label = labeller(variable1 = chemicals_label, variable2 = chemicals_label)) +
  scale_y_continuous(breaks = c(-0.4, 0, 0.4)) + 
  xlab("Chemical") +
  ylab("Externalizing") +
  scale_color_manual(values = c("steelblue3", "salmon3", "springgreen4")) +
  theme_bw(base_family = "AppleGothic") +
  theme(text = element_text(size = 30, color = "grey25"),
      axis.text = element_text(color = "grey25", size = 30),
        strip.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
       legend.position = "bottom")


figure_s3 <- plot_grid(figure_s21,figure_s22, col = 2, align = "h")
figure_s3

Self-Organizing Maps

Figure 4

Self-organizing maps (SOM) cluster star plot reflecting patterns of exposure to per- and polyfluoroalkyl substances (PFAS) and polybrominated diphenyl ethers (PBDEs), measured in serum samples obtained during the first trimester, and association between SOM clusters and internalizing and externalizing behaviors averaged across ages 1-5 years among Black mother-child pairs in metropolitan Atlanta, Georgia.

# Define mixture
df.exp <- mix %>% dplyr::select(PFHXS.log, PFOS.log, PFOA.log, PFNA.log, PBDE47.log, PBDE99.log, PBDE100.log)

df.exp <- df.exp %>% rename(PFHXS = PFHXS.log, PFOS = PFOS.log, PFOA = PFOA.log, PFNA = PFNA.log, "BDE-47" = PBDE47.log, "BDE-99" = PBDE99.log, "BDE-100" = PBDE100.log)

# Explore clustering stats
# grp.str<-grpeval(df.exp, kmx=10)

# Initialize neurons/size of grid
grid <- somgrid(xdim = 3, ydim = 1, topo = "hexagonal")
# Train SOM model
set.seed(123)
som.mix <- som(X = as.matrix(df.exp), grid = grid, rlen = 500)

# Check if number of iterations is ok (should be stable) 
# plot(som.mix, type = "changes")

## Visualize

# N
# plot(som.mix, type = "counts")
# plot(som.mix, type = "mapping", keepMargins = TRUE)

# Codes/weights
par(cex = 8)
plot(som.mix, type = "codes", main = "SOM Clusters")

# Get unit classification
mix$unit_clusters <- som.mix$unit.classif

# get events (weight vector)
# getCodes(som.mix)

# Regression based on unit classifications

# Assign reference group 
mix$unit_clusters <- relevel(as.factor(mix$unit_clusters), ref = "1")

# Internalizing
som_int_reg <- lm(Internalizing ~ unit_clusters + `Mean Child Age`  + `Infant Sex` + `Marital Status` + `Pre-pregnancy BMI` + `Maternal Age` + `Maternal Education` + Parity + `Substance Use`, data = mix)

som_int_reg_df <- tidy(som_int_reg , conf.int = TRUE) %>% filter(term %in% c("unit_clusters2", "unit_clusters3")) %>% dplyr::select(term, estimate, conf.low, conf.high) %>% mutate(outcome = "Internalizing")


# Externalizing 
som_ext_reg <- lm(Externalizing ~ unit_clusters + `Mean Child Age`  + `Infant Sex` + `Marital Status` + `Pre-pregnancy BMI` + `Maternal Age` + `Maternal Education` + Parity + `Substance Use`, data = mix)

som_ext_reg_df <- tidy(som_ext_reg , conf.int = TRUE) %>% filter(term %in% c("unit_clusters2", "unit_clusters3")) %>% dplyr::select(term, estimate, conf.low, conf.high) %>% mutate (outcome = "Externalizing")


som_reg_df <- rbind(som_int_reg_df, som_ext_reg_df )
som_reg_df$outcome <- factor(som_reg_df$outcome, levels = c("Internalizing",  "Externalizing"))

figure_5 <- ggplot(som_reg_df, aes(x = estimate, y = term, color = outcome)) +
  geom_pointrange(aes(xmin = conf.low, xmax = conf.high), size = 4, linewidth = 4) +
  geom_vline(xintercept = 0, color = "grey25", linetype = "dashed", alpha = 0.5) +  
  facet_grid(~outcome, scales = "free") +
  xlab("Beta (95% Confidence Interval)") +
  ylab("Chemical") +
  scale_x_continuous(n.breaks=4) +
  scale_y_discrete(labels = c("Cluster 2", "Cluster 3")) +
  scale_color_manual(values = c("grey25", "grey55")) +
  coord_flip() +
  theme_classic(base_family = "AppleGothic") +
  theme(title = element_text(size = 50, color = "grey25"),
        axis.title.x = element_text(size = 50, color = "grey25"),
        axis.title.y = element_text(size = 50, color = "grey25"),
        axis.text.x = element_text(size = 50, color = "grey25", vjust = 0.6),
        axis.text.y = element_text(size = 50, color = "grey25"),
        strip.background = element_rect(fill="white"),
        strip.text = element_text(colour = 'grey25', face = "bold", size = 50), 
        legend.position = "none", 
        strip.background.x = element_blank(),
        panel.grid.major = element_blank(), panel.grid.minor = element_blank())

figure_5

Table S7

Linear regression associations between self-organizing map (SOM) clusters and internalizing and externalizing behaviors averaged across ages 1-5 years among Black mother-child pairs in metropolitan Atlanta, Georgia.

Models are adjusted for infant sex, mean child age, maternal age, maternal education, marital status, early pregnancy BMI, parity, and substance use.

# Table
table_s8 <- som_reg_df %>%
  mutate_if(is.numeric, round, digits = 2) %>%
  mutate(Beta_CI = paste0(estimate, " (", conf.low, ", ", conf.high, ")"))  %>%
 mutate(term = recode(term, 
                     "unit_clusters2" = "Cluster 2",
                     "unit_clusters3" = "Cluster 3")) %>%
  dplyr::select(term, outcome, Beta_CI) %>%
  pivot_wider(names_from = outcome, values_from = (Beta_CI)) %>%
  flextable() %>%
  add_body_row(values = c("" ,"Beta (95% CI)", "Beta (95% CI)")) %>%
  font(fontname = "Tahoma") %>%
  set_header_labels(term = "", sex = "Sex") %>%
  autofit() %>%
  theme_alafoli()

table_s8

Internalizing

Externalizing

Beta (95% CI)

Beta (95% CI)

Cluster 2

0.41 (-0.01, 0.83)

0.5 (0.09, 0.91)

Cluster 3

0.44 (0.08, 0.81)

0.54 (0.18, 0.89)